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Effects  of  stencil  width  on  surface  ocean  geostrophic  velocity 
and  vorticity  estimation  from  gridded  satellite  altimeter  data 

Brian  K.  Arbic,1 * *  Robert  B.  Scott,2,3  Dudley  B.  Chelton,4  James  G.  Richman,5 
and  Jay  F.  Shriver5 
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[i]  This  paper  examines  the  effect  of  “stencil  width”  on  surface  ocean  geostrophic  velocity 
and  vorticity  estimated  from  differentiating  gridded  satellite  altimeter  sea  surface  height 
products.  In  oceanographic  applications,  the  value  of  the  first  derivative  at  a  central 
grid  point  is  generally  obtained  by  differencing  the  sea  surface  heights  at  adjacent  grid 
points.  This  is  called  a  “three-point  stencil  centered  difference”.  Here  the  stencil  width 
is  increased  from  three  to  five,  seven,  and  nine  points,  using  well-known  formulae  from 
the  numerical  analysis  literature.  The  discrepancies  between  velocities  computed  with 
successive  stencils  decreases  with  increasing  stencil  width,  suggesting  that  wide  stencil 
results  are  more  reliable.  Significant  speed-dependent  biases  (up  to  10-20%)  are 
found  between  results  computed  from  three-point  stencils  versus  those  computed  from 
wider  stencils.  The  geostrophic  velocity,  and  the  variance  of  geostrophic  velocity,  are 
underestimated  with  thin  stencils.  Similar  results  are  seen  in  geostrophic  velocities 
computed  from  high-resolution  model  output.  In  contrast  to  the  case  when  three-point 
stencils  are  used,  wider  stencils  yield  estimates  of  the  anisotropy  of  velocity  variance  that  are 
insensitive  to  the  differences  in  grid  spacing  between  two  widely  used  altimeter  products. 

Three-point  stencils  yield  incorrect  anisotropies  on  the  1/4°  anisotropic  AVISO  grid; 
we  recommend  the  use  of  7-point  stencils.  Despite  the  demonstrated  inadequacies  of  the 
three-point  stencils,  the  conclusions  of  earlier  studies  based  on  them,  that  the  zonally 
averaged  midlatitude  eddy  kinetic  energy  field  is  nearly  isotropic,  are  found  to  pertain 
also  with  wider  stencils.  Finally,  the  paper  also  examines  the  strengths  and  limitations 
of  applying  noise-suppressing  differentiators,  versus  classic  centered  differences, 
to  altimeter  data. 

Citation:  Arbic,  B.  K.,  R.  B.  Scott,  D.  B.  Chelton,  J.  G.  Richman,  and  J.  F.  Shriver  (2012),  Effects  of  stencil  width  on  surface 
ocean  geostrophic  velocity  and  vorticity  estimation  from  gridded  satellite  altimeter  data,  J.  Geophys.  Res.,  117 ,  C03029, 
doi:  10.1 029/20 11JC007367. 


1.  Introduction 

[2]  Satellite  altimeter  measurements  of  sea  surface  height 
[Fu  and  Cazenave ,  2001]  have  revolutionized  physical 
oceanography,  enabling  and  enhancing  numerous  avenues 
of  study.  To  name  just  three  of  many  applications,  satellite 
altimetry  data  has  been  utilized  to  estimate  sea  surface  height 
and  kinetic  energy  spectra  [Stammer,  1997;  Xu  and  Fu , 
2011],  reveal  the  inverse  kinetic  energy  cascade  of  surface 


’Department  of  Earth  and  Environmental  Sciences,  University  of 
Michigan,  Ann  Arbor,  Michigan,  USA. 

institute  for  Geophysics,  Jackson  School  of  Geosciences,  University  of 
Texas  at  Austin,  Austin,  Texas,  USA. 

^D6partcmcnt  dc  Physique  and  LPO,  University  dc  Bretagne 
Occidental,  Brest,  France. 

4CoIIcgc  of  Oceanic  and  Atmospheric  Sciences,  Oregon  State 
University,  Corvallis,  Oregon,  USA. 

5Oceanography  Division,  Naval  Research  Laboratory,  Stcnnis  Space 
Center,  Mississippi,  USA. 

Copyright  20 1 2  by  the  American  Geophysical  Union. 

0148-0227/12/201 1JC007367 


oceanic  flows  [Scott  and  Wang ,  2005],  and  document  the 
propagation  paths  of  mesoscale  eddies  [Chelton  et  al ,  2007, 
201 1].  A  recent  brief  review  of  applications  of  altimeter  data 
is  given  by  Scott  et  al  [2010]. 

[3]  Two  important  products  of  satellite  altimeter  data  are 
the  ocean  surface  velocity  and  vorticity  fields,  extracted  via 
the  geostrophic  relations.  Geostrophic  velocities  are  obtained 
from  first  derivatives  of  the  sea  surface  height  77: 


u 


g&n 

fdy' 


(1) 


V 


_gfy 

fdx' 


(2) 


where  u  is  zonal  velocity  (velocity  in  the  east-west  direction, 
with  eastward  velocities  being  positive),  v  is  meridional 
velocity  (velocity  in  the  north-south  direction,  with  north¬ 
ward  velocities  being  positive),  g  is  gravitational  accelera¬ 
tion,  /  is  the  Coriolis  parameter,  x  is  the  zonal  spatial 
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coordinate,  and  y  is  the  meridional  spatial  coordinate.  Geo- 
strophic  vorticity  (  is  obtained  from  first  derivatives  of  the 
velocity  field,  and  therefore  second  derivatives  of  sea  surface 
height: 

dv_5u_£/dV  aVv  _g/?dr; 

4 ~dx  dy~fWYdf)  p  dy'  K) 

where  (3  =  d  fldy  is  the  planetary  vorticity  gradient.  The  ratio 
of  the  last  term  on  the  right-hand  side  in  (3)  to  the  second 
derivative  terms  scales  as  (3L!f9  where  I  is  a  typical  horizontal 
length  scale.  Under  the  traditional  “beta  plane  approxima¬ 
tion”,  this  ratio  is  small  in  midlatitudes  [e.g.,  Pedlosky ,  1987; 
Vallis,  2006]. 

[4]  Here  we  examine  the  effects  of  “stencil  width”  on  the 
estimation  of  geostrophic  velocities  and  vorticities  from 
gridded  altimeter  sea  surface  height  products.  By  “stencil 
width”  we  mean  the  number  of  grid  points  utilized  to  estimate 
the  finite  difference  approximation  to  the  derivative  on  a  grid. 
As  we  will  discuss,  varying  the  stencil  width  is  another  way 
of  varying  the  accuracy  of  the  derivative  estimate.  Our  focus 
will  be  on  the  velocity  fields,  but  we  will  include  a  brief 
discussion  of  the  vorticity  field.  In  the  vast  majority  of 
oceanographic  applications,  geostrophic  velocities  are  esti¬ 
mated  from  a  “three-point  stencil  centered  difference”,  in 
which  sea  surface  height  values  at  adjacent  grid  points  are 
differenced  to  determine  velocities  at  a  central  grid  point. 
However,  inspection  of  the  substantial  literature  on  numeri¬ 
cal  methods  [e.g.,  Strikwerda ,  2004;  Mathews  and  Fink , 
2004]  reveals  well-known  formulae  for  computing  first 
derivatives  of  fields  via  centered  differences  taken  with  wider 
stencils,  for  instance,  5-point,  7-point,  and  9-point  stencils. 
(We  also  found  the  following  websites  helpful:  http:// 
www.holoborodko.com/pavel/numerical_methods/numerical- 
derivative/,  which  was  especially  helpful  and  is  herein¬ 
after  referred  to  as  Holoborodko  (201 1);  http ://en. wikipedia. 
org/wiki/Five-point_stencil;  http://reference.wolfram.com/ 
mathematica/tutorial/NDSolvePDE.html#c:4.)  It  is  of  inter¬ 
est  to  determine  whether  usage  of  wider  stencils  would  alter 
the  velocity  and  vorticity  estimates.  As  we  shall  see  from 
the  numerical  analysis  literature,  derivatives  estimated  with 
wider  stencils  feature  reduced  errors  over  those  computed 
with  3-point  stencils. 

[5]  Since  the  errors  in  derivative  estimates  are  a  function 
of  grid  spacing,  it  is  also  of  interest  to  determine  how  stencil 
width  affects  velocity  estimates  computed  on  grids  of  dif¬ 
fering  grid  spacings.  For  example,  we  will  examine  deriva¬ 
tives  computed  on  the  widely-used  isotropic  1/3°  Mercator 
grid  reference  product  distributed  by  the  Archiving,  Valida¬ 
tion,  and  Interpretation  of  Satellite  Oceanographic  (AVISO) 
gridded  altimeter  data  [Le  Traon  et  al9  1998;  Ducet  et  al. , 
2000],  versus  derivatives  computed  on  the  (also  widely  used) 
anisotropic  1/4°  latitude-longitude  grid  reference  product  put 
out  by  AVISO.  By  “isotropic  grid”,  we  mean  a  grid  in  which 
the  spacing  in  kilometers  is  the  same  in  the  east-west  and 
north-south  directions.  This  is  the  case  for  the  1/3°  Mercator 
AVISO  grid,  in  which  the  zonal  grid  spacing  is  fixed  in 
degrees  of  longitude  and  the  meridional  grid  spacing 
decreases  in  degrees  of  latitude  with  the  cosine  of  latitude 
so  as  to  match  the  zonal  grid  spacing  in  kilometers.  On 
an  “anisotropic  grid”,  such  as  the  1/4°  AVISO  latitude- 
longitude  grid,  the  grid  spacing  in  the  north-south  direction 


is  the  same  across  the  entire  grid,  while  grid  spacing  in  the 
east-west  direction  is  fixed  in  degrees  of  longitude  but 
decreases  in  kilometers  with  the  cosine  of  latitude.  Since  the 
1/4°  product  is  simply  a  bi-linear  interpolation  of  the  1/3° 
product,  we  might  expect  that  derivatives  calculated  on  the 
two  grids  should  yield  very  similar  results.  We  will  show  that 
this  is  the  case  only  if  wide  stencils  are  utilized.  Since  the 
derivative  in  the  zonal  direction  yields  meridional  velocity 
while  the  derivative  in  the  meridional  direction  yields  zonal 
velocity,  this  discussion  will  be  useful  for  the  estimation  of 
anisotropy  in  the  oceanic  kinetic  energy  field  [Ducet  et  al ., 
2000;  Scott  et  al .,  2008;  Scharffenberg  and  Stammer ,  2010]. 

[6]  We  will  also  test  whether  the  behaviors  seen  when 
differentiating  sea  surface  heights  in  altimetric  data  sets  are 
seen  as  well  when  differentiating  sea  surface  heights  in  high- 
resolution  numerical  models.  Models  and  altimetric  data  sets 
exhibit  different  types  of  errors,  hence  a  consistency  in 
results  obtained  from  the  two  sources  would  demonstrate  that 
the  impact  of  stencil  width  is  not  simply  an  artifact  of  the 
particular  nature  of  gridded  altimeter  products.  Here  we  uti¬ 
lize  results  from  NLOM,  the  Naval  Research  Laboratory 
Layered  Ocean  Model  [Hurlburt  and  Thompson ,  1980; 
Wallcraft  et  al.9  2003].  NLOM  is  run  as  a  data-assimilative 
nowcast/forecast  model  by  the  United  States  Navy  [Shriver 
et  al .,  2007].  Here  for  simplicity  we  utilize  a  snapshot  from 
a  forward  (non-assimilative)  run  of  NLOM.  Since  the  model 
is  on  a  much  higher  horizontal  resolution  grid  than  is  the 
altimeter  data,  we  can  use  either  subsampled  or  smoothed 
versions  of  the  model  to  further  test  the  effect  of  stencil 
widths  for  differing  horizontal  grid  resolutions. 

[7]  Most  of  the  paper  is  about  “classic”  centered  differ¬ 
ences,  which  operate  well  on  noiseless  data  or  on  numerical 
models.  Later  in  the  paper  we  will  examine  the  impact  of  data 
noise  on  velocity  estimates.  We  will  also  discuss  application 
of  noise-suppressing  differentiators,  which  have  been  widely 
discussed  in  the  chemical  and  signal  processing  literature 
[Savitzky  and  Golay ,  1964;  Steiner  et  al.9  1972;  Gorry ,  1990; 
Luo  et  al.9  2005;  Holoborodko,  2011],  to  altimeter  data. 

[8]  This  paper  is  organized  as  follows.  In  section  2,  we  list 
the  formulae  for  first  and  second  derivatives  computed  as  3-, 
5-,  7-,  and  9-point  stencil  “classic”  centered  differences, 
taken  from  standard  references  in  the  numerical  finite  dif¬ 
ference  literature.  We  also  list  the  leading  order  error  terms 
associated  with  the  derivative  estimates.  In  section  3,  we 
show  the  results  of  a  “von  Neumann  analysis”  [ LeVeque , 
2007],  which  demonstrates  that  the  deviations  of  classic 
centered-difference  estimates  of  derivatives  from  estimates 
made  by  an  “ideal  differentiator”  (to  be  explained  in 
section  3)  decrease  as  stencil  width  increases.  We  utilize 
results  from  an  idealized  quasi-geostrophic  turbulence  model 
in  this  exercise.  Since  these  deviations  depend  on  length 
scale,  we  include  and  discuss  a  wave  number  spectrum  of  the 
AVISO  data  in  section  3.  We  then  show,  in  section  4,  that  the 
differences  between  geostrophic  velocities  computed  with 
successively  wider  stencils  decreases  with  increasing  stencil 
width.  Section  4  utilizes  fields  of  sea  surface  height  taken 
from  the  weekly  1/3°  and  1/4°  AVISO  products  as  well  as 
from  NLOM.  Section  5  presents  a  brief  discussion  of  the 
alterations  to  the  classic  centered  difference  formulae  when 
non-uniform  grid  spacing  is  taken  into  account  This  is  a 
small  but  measurable  effect,  relevant  for  zonal  velocities 
computed  on  the  1/3°  AVISO  grid.  In  section  6,  we  present  a 
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Table  1.  Formulae  for  Classic  Centered  Differences 

N *  Formulae  Eb 


3 

5 

7 

9 

3 

5 

7 
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First  Derivative c 
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“Value  of  A,  the  width  of  the  stencil, 
heading-order  error  term  E  in  estimate. 

cFormulac  for  Y-point  classic  centered  difference  estimate  of  first  derivative  of  sea-surface  elevation  rj. 
dFormulac  for  A-point  classic  centered  difference  estimate  of  second  derivative  of  sea-surface  elevation  r /. 


brief  discussion  of  geostrophic  vorticity  computed  from 
AVISO  data  with  classic  centered  differences.  This  dis¬ 
cussion  continues  to  demonstrate  the  advantages  of  wider 
stencils.  In  section  7,  we  present  a  stochastic  model  which 
qualitatively  explains  some  of  the  “speed-dependent  biases” 
seen  in  the  three-point  stencil  results  of  section  4.  In 
section  8,  we  show  that  consistent  estimates  of  the  anistropy 
of  kinetic  energy  variance  on  anisotropically  versus  isotrop¬ 
ically  gridded  altimeter  products  can  be  obtained  if  and  only 
if  wide  stencils  are  utilized.  Section  9  examines  the  impact  of 
noise  on  the  derivative  estimates,  and  introduces  some  noise- 
suppressing  differentiators  taken  from  Holoborodko  (2011). 
We  argue  that  for  current-generation  altimeter  data  (AVISO), 
the  merits  of  suppressing  noise  at  small  scales  are  out¬ 
weighed  by  the  need  to  retain  near  “ideal  differentiation”  at 
larger  scales,  which  the  wider  stencil  classic  centered  differ¬ 
ences  provide.  In  section  10  we  present  a  brief  summary  and 
discussion  of  our  results. 

2.  Formulae  for  N-Point  Stencil  Classic  Centered 
Differences 

[9]  For  geostrophic  velocity  and  vorticity  estimation,  we 
need  to  compute  partial  first  and  second  derivatives  of  the  sea 
surface  height  7]  with  respect  to  the  east-west  coordinate  x 
and  north-south  coordinate  y.  Let  r  be  a  generic  coordinate, 
representing  either  x  or  y.  Suppose  further  that  we  take  the 
grid  spacing  of  an  altimeter  product  to  be  h.  For  the  moment 
we  assume  uniform  grid  spacing.  The  3-,  5-,  7-,  and  9-point 
stencil  estimates  of  the  first  derivative  utilize  1,  2,  3,  and 
4  grid  points,  respectively,  on  either  side  of  the  grid  point  at 
which  the  derivative  is  required.  For  instance,  the  3-point 
stencil  utilizes  the  value  of  7]  at  r  +  h  and  r-h  to  estimate  the 
first  derivative  at  r,  the  5-point  stencil  utilizes  the  value  of  77 
at  r  +  2h,  r  +  h,  r-h,  and  r-2h  to  estimate  this  derivative,  and 
so  on.  Table  1  gives  the  formulae  for  the  first  and  second 
partial  derivatives  of  7]  computed  as  classic  centered  differ¬ 
ences  on  the  3-,  5-,  7-,  and  9-point  stencils,  as  well  as  the 
leading-order  error  term  E  in  the  estimate  [e.g.,  Strikwerda , 
2004;  Mathews  and  Fink ,  2004]  (http://www.holoborodko. 
com/pavel/numerical  methods/numerical-derivative/,  http://en. 
wikipedia.org/wiki/Five-point_stencil,  http://reference. wolfram. 


com/mathematica/tutorial/NDSolvePDE.html#c:4).  Note  that 
superscripts  in  parentheses  denote  order  of  derivative.  Thus, 
for  example,  Ty\r)  denotes  the  ninth  derivative  of  r£r).  Note 
also  that  the  denominator  in  the  error  term  increases 
with  stencil  width.  Therefore  we  expect  the  differences  in 
estimates  made  with  successive  stencils  to  decrease  as  the 
stencil  widens. 

[10]  For  the  benefit  of  readers  not  familiar  with  the 
numerical  finite  difference  methods  literature,  the  derivation 
of  the  first  derivative  and  error  formulae  is  sketched  out 
below.  (For  more  detail,  the  reader  can  consult  the  numerical 
analysis  references  listed  above,  or  Cushman-Roisin  and 
Beckers  [2010,  chapter  1].)  Using  the  9-point  stencil  com¬ 
putation  as  an  example,  we  write  the  derivative  of  sea  surface 
height  at  a  given  grid  point  as  a  linear  combination  of  sea 
surface  heights  at  the  four  nearest  grid  points  on  each  side, 

^~^Y^ciV(r  +  nh),  n  =  -4,  -3,  -2,-1, 1,2, 3, 4.  (4) 
For  each  value  of  n ,  we  write  a  Taylor  series 

.  WV3)M  ,  WW) 

3J  91  ’  w 

where  we  truncate  the  series  at  the  ninth-derivative  term. 
Inserting  the  Taylor  series  into  (4)  yields 

+  f|  (c„)r/<l>(r)  +  F2(c„)hr/<2>(r)  +  ... 

+  F9(c„)hsrjm(r),  (6) 

where  F0(c„)f  Fi(c„),  ...  F9(cn)  are  simple  linear  combi¬ 
nations  of  the  coefficients  cn.  In  order  that  only  the  first 
derivative  is  retained  in  the  right-hand  side  of  (6),  we  want 
F\(cn)  to  be  unity,  and  Fn  for  n  £  1  to  be  zero  for  as  many  n 
as  possible.  For  a  9-point  stencil  we  are  able  to  set  the  Fn  to 
zero  for  n  £  1 ,  all  the  way  up  to  F8.  Setting  F0(c„)  =  F2(cn)  = 
F4(cn)  =  F6(c„)  =  0  yields  c_„  =  —  cn  for  n  ~  1,  2,  3,  4.  This 
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(cycles  per  degree  longitude)  (cycles  per  degree  longitude) 


Figure  1.  The  zonal  wave  number  spectra  of  sea  surface  height  (left)  from  the  AVISO  1/3°  Mercator 
gridded  product  and  (right)  from  the  AVISO  1/4°  product.  See  text  for  discussion  of  normalization.  The 
vertical  dashed  and  dotted  lines  correspond  to  wavelengths  of  3°  and  2°,  respectively. 


implies  in  turn  that  F%(cn)  also  equals  zero.  We  are  then  left 
with 

^  ~  F,  (c«)*7(,)M  +  F3(c„)h2rj{3)(r)  +  F5(cn)h4rf5)(r) 

+  F7(cn)h6rj^(r)  +  F9(cn)h*VW(r)y  (7) 

where  now  each  Fn  depends  only  on  the  four  constants  C\,  C2 , 
C3,  and  c4,  since  we  have  written  in  terms  of  cn  for  n  =  1, 

2,  3,  4.  Setting  Fx(cn)  =  1  and  F3(c„)  =  F5(c„)  =  F7(c„)  =  0 
provides  four  equations  which  can  be  solved  for  the  four 
constants  cm  n  =  1,  2,  3,  4  as  shown  in  Table  1.  We  are  then 
left  with  F9(cn)hsrf9\r)  as  the  leading-order  error  term.  The 
formulae  for  second  derivatives  are  obtained  in  like  manner, 
the  main  difference  being  that  the  value  of  rj  at  the  central  grid 
point  r  is  utilized: 

| ^E^  +  nh),  "  =  -4,-3, -2, -1,0, 1,2, 3, 4.  (8) 

3.  The  von  Neumann  Analysis  and  Spectral 
Content  of  AVISO  Data 

[11]  We  can  anticipate  some  of  the  results  of  using  wider 
stencils  by  following  the  “von  Neumann  analysis”  given  in 
section  9.6  of  LeVeque  [2007].  This  analysis  allows  us  to 
examine  how  accurate  the  stencil  approximations  for  deri¬ 
vatives  are,  as  a  function  of  the  differing  length  scales  (wave 
numbers)  contained  in  the  spatially  varying  sea  surface 
height,  rj(r).  We  focus  on  the  first  derivative  in  this  section. 

[12]  Because  of  the  importance  of  the  wave  number  con¬ 
tent  of  the  AVISO  signal  to  the  von  Neumann  analysis,  we 
display  in  Figure  1  the  wave  number  spectrum  of  AVISO. 
The  figure  displays  the  zonal  wave  number  spectrum  of  sea 
surface  height  for  wave  number  in  units  of  cycles  per  degree 
of  longitude,  computed  along  1 5  zonal  sections  in  the  North 
and  South  Pacific  [see  Chelton  et  al ,  2011,  Figure  A2]. 
Spectra  from  both  the  1/3°  Mercator  grid  AVISO  product, 
and  the  1/4°  AVISO  product,  are  shown.  The  spectra  have 


been  normalized  so  that  each  spectrum  has  the  same  variance 
integrated  over  wavelengths  shorter  than  3°  (wave  numbers 
higher  than  0.333  cycles  per  degree  of  longitude).  As  dis¬ 
cussed  by  Chelton  et  al  [2011],  while  spectra  for  wave¬ 
lengths  larger  than  2°-3°  display  substantial  variance  from 
one  section  to  the  next  (note  the  scatter  of  the  different  curves 
in  the  left-hand  side  of  both  Figures  1  (left)  and  1  (right)), 
spectra  for  the  different  sections  over  wavelengths  shorter 
than  2°-3°  lie  almost  on  top  of  each  other.  Chelton  et  al 
[2011]  argue  that  this  demonstrates  that  the  filtering  inher¬ 
ent  in  the  creation  of  AVISO  products  removes  most  of  the 
variance  in  wavelengths  shorter  than  2°-3°.  We  will  refer 
to  the  2°-3°  scale  as  the  wavelength  resolution  limit.  The 
bi-linear  interpolation  introduces  a  sidelobe  with  peak  at 
1 .333  cycles  per  degree  into  the  zonal  wave  number  spectrum 
of  kinetic  energy  on  the  1/4°  grid  (compare  Figure  1  (right)  to 
Figure  1  (left)).  The  sidelobe  is  absent  in  spectra  computed 
on  the  1/3°  grid  [Chelton  et  al ,  2011]  and  is  therefore  an 
artifact  introduced  by  interpolation. 

[13]  Since  the  sea  surface  height  field  rj  on  a  grid  r  can  be 
written  as  a  Fourier  sum 

k 

where  /  =  ,  it  suffices  to  examine  the  derivative  of 

the  function  elkr.  This  function  is  an  eigenfunction  of  the 
“ideal  differentiator”  with  eigenvalue  A  =  ik,  since  ^  el*r  = 
ike lkr.  Let  us  now  examine  the  stencil  approximation  for  the 
derivative  of  77  =  elk{r  +  nA),  the  discretized  equivalent  of  elkr. 
As  before,  n  is  the  index  of  the  grid  point  number  relative 
to  the  central  grid  point  r,  and  h  is  the  grid  point  spacing. 
The  wave  number  k  =  where  A  is  wavelength.  Let  A  be 
the  2m  +  1  point  stencil  centered  difference  operator.  Then 
we  have 

&n(r)  =  j  5Z  c»rAr+nh)  (n^°).  0°) 
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or,  since  c_n  =  —  cn , 


Figure  2.  Response  of  an  ideal  differentiator  and  of  3-,  5-, 
7-,  and  9-point  stencils,  as  a  function  of  wave  number,  on 
(a)  1/3°  AVISO  grid  and  (b)  NLOM  grid.  The  grid  spacings 
used  for  Figures  2a  and  2b  are  37  and  4.9  km,  respectively, 
and  the  response  is  the  absolute  value  \A\  of  the  eigenvalue 
of  the  von  Neumann  analysis.  See  text  for  more  details. 
Vertical  lines  are  drawn  at  the  wave  number  (40  km)  \ 
corresponding  approximately  to  the  2-3°  wavelength  resolu¬ 
tion  limitation  of  the  AVISO  sea  surface  height  fields. 


At ?(r)  =  !£>, -  «'*('—*> ] .  (11) 

Therefore 

Ar?(r)  =  X-jrk  c„[2isin(nA*)].  (12) 

rt=*l 

Therefore  eirk  is  an  “eigengrid function”  of  the  operator  A, 
with  eigenvalue 

A  =  c„sin(nhk).  (13) 

»»i 

For  example,  for  a  3-point  stencil,  m  =  1,  C\  —  and  A  = 
jj  sin (hk).  For  small  hk  we  can  write  A~  j [  ^hk  —  |  (hk)*  + . . .J  ~ 

ik j^l  —  \(hk)2  -f  ...j.  This  agrees  with  the  eigenvalue  ik  of 
J:  to  order  (hk)2. 

[14]  In  Figure  2  we  plot  the  results  of  the  von  Neumann 
analysis  for  classic  centered  differences.  We  plot  the 
“response”  IAI  for  the  3-point,  5-point,  7-point,  and  9-point 
stencils,  as  a  function  of  wave  number  k.  We  also  plot  k ,  the 
value  of  IAI  for  the  ideal  differentiator  The  “response”, 
represented  by  the  vertical  axis  in  figures  like  this,  is  some¬ 
times  referred  to  as  the  “filter  transfer  function”  in  the  signal 
processing  literature.  In  Figure  2a  we  set  h  -  37  km,  the 
equatorial  grid  spacing  in  the  x-  and  y-directions  on  the  1/3° 
AVISO  grid.  In  Figure  2b  we  set  h  =  4.9  km,  the  equatorial 
grid  spacing  in  the  x-direction  on  the  NLOM  grid.  As  a  result, 
note  that  the  x-axes  in  Figures  2a  and  2b  differ  from  each 
other.  In  both  plots  we  draw  in  (40  km)-1,  the  approximate 
wave  number  corresponding  to  the  ~2°  -  3°  wavelength 
resolution  limitation  of  the  AVISO  sea  surface  height  fields 
(Figure  1)  [Chelton  et  al ,  201 1],  as  a  vertical  line.  Figure  2a 
shows  that  for  wave  numbers  just  below  (length  scales  just 
above)  the  resolution  limitation,  the  3-point  stencil  deviates 
significantly  from  the  ideal  response.  This  implies  that  an 
important  part  of  the  signal  is  lost  when  3-point  stencils  are 
used.  The  wider  stencils  do  not  deviate  from  the  ideal  dif¬ 
ferentiator  until  larger  wave  numbers  (smaller  scales)  are 
reached.  From  this  we  anticipate  that  there  will  be  substantial 
differences  between  derivatives  computed  using  3-point 
versus  wider  stencils  on  the  1/3°  AVISO  grid  (and  on  the  1/4° 
AVISO  grid,  since  the  grid  spacings  are  comparable).  In 
Figure  2b,  because  of  the  smaller  grid  spacing,  the  wave 
number  (40  km)-1  lies  far  below  the  wave  numbers  where 
the  3-point  stencil  and  other  stencils  deviate  from  the  ideal 
differentiator.  For  features  with  length  scales  larger  than 
about  40  km,  we  therefore  expect  stencil  width  to  make  less 
of  a  difference  on  the  high-resolution  NLOM  grid  than  on  the 
coarser  AVISO  grid. 

[15]  To  illustrate  the  implications  of  the  von  Neumann 
analysis,  we  utilized  results  from  the  idealized  two-layer 
doubly  periodic  pseudo-spectral  quasi-geostrophic  turbulence 
simulations  of  Arbic  and  Flierl  [2004].  Here  we  used  the 
snapshot  in  their  Figure  8d.  The  advantage  of  using  the 
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Figure  3.  Differences  between  zonal  velocity  component  u  computed  using  3-point,  5-point,  7-point,  and 
9-point  stencils,  and  u  computed  using  an  “ideal  differentiator”  spectral  technique  (see  text).  Velocities 
are  computed  from  a  snapshot  of  the  upper  layer  stream  function  in  an  idealized  doubly  periodic  quasi- 
geostrophic  turbulence  model.  The  y-axis  of  the  scatterplot  is  the  difference  between  u  computed  on 
stencils  and  u  computed  from  the  ideal  differentiator,  while  the  x-axis  represents  u  computed  from  the  ideal 
differentiator.  Red,  green,  black,  and  magenta  dots  represent  differences  computed  using  3-point,  5-point, 
7-point,  and  9-point  stencils,  respectively.  Both  axes  are  normalized  by  the  imposed  time-mean  flow  in  the 
idealized  model,  which  was  taken  to  be  1  cms"1  in  the  work  of  Arbic  and  Flierl  [2004]. 


idealized  model  in  this  illustration  is  that  velocities  are 
computed  in  this  model  via  the  “ideal  differentiator”  of  the 
von  Neumann  analysis.  (Indeed,  this  property  is  one  reason 
pseudo-spectral  models  have  enjoyed  a  long  usage.)  In 
pseudo-spectral  models,  the  model  stream  function  is  Fourier 
transformed,  multiplied  by  -il  or  ik,  and  then  inverse  trans¬ 
formed  to  obtain  u  and  v,  respectively,  where  /  is  the  merid¬ 
ional  wave  number  and  k  is  the  zonal  wave  number.  Since 
the  idealized  model  is  doubly  periodic,  no  tapering  near  the 
edges  is  required  to  perform  this  computation,  in  contrast 
to  spectral  computations  done  in  more  realistic  domains. 

[16]  Figure  3  displays  scatterplots  of  the  differences 
ustencii  —  uideal  between  zonal  velocity  ustencil  computed  from 
the  idealized  model  using  3-point,  5-point,  7-point,  and 
9-point  stencils,  and  uideai  computed  using  the  “ideal  differ¬ 
entiator”  spectral  technique  actually  employed  in  the  ideal¬ 
ized  model.  The  differences  ustencii  —  uideai  are  plotted  on  the 
y-axis  of  Figure  3,  while  the  x-axis  displays  uideaJ  computed 
from  the  ideal  differentiator.  The  differences  decrease  as 
stencil  width  increases,  consistent  with  the  von  Neumann 
analysis  discussed  above.  Because  the  dynamics  in  this  par¬ 
ticular  simulation  of  the  idealized  model  are  nearly  isotropic, 
and  because  the  grid  is  also  isotropic,  results  for  the  meridi¬ 
onal  velocity  v  are  nearly  identical  to  those  in  Figure  3, 
and  are  not  shown  for  the  sake  of  brevity.  Because  veloc¬ 
ities  computed  from  the  9-point  stencil  lie  closer  to  velocities 
computed  from  the  ideal  differentiator  than  do  velocities 
computed  from  narrower  stencils,  we  will  take  the  9-point 
stencil  results  as  our  standard  in  what  follows.  Ideally,  we 


will  expect  7-point  stencil  results  to  lie  closer  to  9-point 
stencil  results  than  5-point  stencil  results  do,  and  5-point 
stencil  results  to  lie  closer  to  9-point  stencil  results  than 
3-point  stencil  results  do. 

4.  Effect  of  Stencil  Width  on  Geostrophic  Velocity 
Estimates  From  Altimeter  Data  and  Realistic 
Models 

[n]  In  this  section  we  examine  the  impact  of  stencil  width 
on  geostrophic  velocity  estimates  computed  using  classic 
centered  differences  from  a  single  snapshot  of  sea  surface 
height.  We  utilize  the  AVISO  weekly  reference  (two  satel¬ 
lite)  product  for  November  5,  2008.  Both  zonal  and  meridi¬ 
onal  velocity  components  are  computed,  and  results  from 
differentiating  on  both  the  1/3°  and  1/4°  AVISO  grids  are 
shown.  The  1/3°  AVISO  grid  has  915  by  1080  grid  points  in 
the  north-south  and  east-west  directions,  respectively,  while 
the  1/4°  AVISO  grid  has  721  by  1440  grid  points.  We  also 
utilize  the  February  15,  2002  snapshot  of  sea  surface  height 
from  NLOM,  which  has  4384  by  8192  grid  points  on  a  near- 
global  domain  (72°S  to  65°N).  The  NLOM  grid  point  spac¬ 
ing  is  1/32°  in  the  north-south  direction  and  45/1024°  in  the 
east-west  direction.  In  this  and  succeeding  results  sections, 
we  choose  to  emphasize  the  deep  ocean,  in  both  the  AVISO 
and  NLOM  calculations.  To  do  so,  we  display  only  velocities 
computed  at  grid  points  for  which  the  seafloor  depth  exceeds 
1000  m.  Furthermore,  where  comparisons  are  made  between 
u  and  v  or  between  velocity  components  computed  from 
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Figure  4.  Scatterplots  of  velocity  components  computed  from  3-point,  5-point,  and  7-point  stencils  minus 
components  computed  from  9-point  stencils  (y-axis),  plotted  versus  components  computed  from  9-point 
stencils  (x-axis).  Red,  green,  and  black  dots  denote  results  with  3-point,  5-point,  and  7-point  stencils, 
respectively,  (a)  Zonal  velocity  u  computed  on  the  1/3°  Mercator  grid  reference  AVISO  product  for 
November  5, 2008.  (b)  As  in  Figure  4a  but  for  meridional  velocity  v.  (c)  As  in  Figure  4a  but  for  zonal  veloc¬ 
ity  u  computed  on  the  1/4°  latitude-longitude  grid  reference  AVISO  product  for  November  5, 2008.  (d)  As 
in  Figure  4c  but  for  meridional  velocity  v. 


different  stencil  widths,  we  only  display  points  for  which  all 
relevant  velocity  components  are  defined.  For  wider  stencils, 
we  will  lose  points  near  boundaries,  where  a  full  span  of  the 
stencil  does  not  exist  because  of  land. 

[is]  Figure  4  shows  scatterplots  of  geostrophic  velocity 
components  computed  from  the  aforementioned  AVISO 
products.  In  like  manner  to  Figure  3,  we  plot  unarrow  —  u9pt 
versus  u9pt  and  vnarrow  —  v9pl  versus  v9ph  where  the  “9/?f” 
subscript  denotes  a  9-point  stencil  difference  and  “ narrow ” 
denotes  either  a  3-point,  5-point,  or  7-point  stencil  difference. 
As  anticipated  in  the  previous  section,  whereas  the  differ¬ 
ences  between  7-point  and  9-point  stencil  results  (black  dots) 
lie  close  to  zero  for  all  grid  points,  the  differences  between 
3-point  and  9-point  stencil  results  (red  dots)  lie  farther  from 
zero  (display  much  more  scatter).  The  scatter  is  not  distrib¬ 
uted  around  zero  evenly,  but  rather  reveals  what  we  will  call 


a  “speed-dependent  bias”.  By  this  we  mean  that  where  the 
9-point  stencil  results  are  negative,  the  3-point  minus  9-point 
stencil  results  are  positive,  and  vice  versa.  The  difference 
between  5-point  and  9-point  stencil  results  (green  dots) 
displays  less  speed-dependent  bias  than  the  3 -point  minus 
9-point  results  but  more  speed-dependent  bias  than  the 
7-point  minus  9-point  results.  The  speed-dependent  bias  at 
individual  grid  points  can  be  as  large  as  10%  or  more.  These 
behaviors  are  seen  for  both  the  u  and  v  components  of 
velocity,  and  for  both  velocities  computed  on  the  1/3° 
AVISO  grid  and  the  1/4°  AVISO  grid.  A  stochastic  model  for 
the  speed-dependent  bias  will  be  presented  in  section  7. 

[19]  Figure  5  demonstrates  that  similar  behaviors  are  seen 
in  geostrophic  velocities  computed  through  differentiation  of 
the  sea  surface  heights  in  NLOM.  Here  we  compute  geo¬ 
strophic  velocities  on  both  the  native  4384  by  8192  NLOM 
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(a)  u,  native  NLOM  grid  (b)  v,  native  NLOM  grid 


(c)  u,  NLOM  grid  decimated  by  8  (d)  v,  NLOM  grid  decimated  by  8 
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Figure  5.  As  in  Figure  4,  but  using  February  15,  2002  output  from  the  Naval  Research  Laboratory 
Layered  Ocean  Model  (NLOM)  rather  than  gridded  altimeter  data,  (a)  Zonal  velocity  u  and  (b)  meridional 
velocity  v  computed  on  NLOM  output  on  the  native  4384  by  8192  grid,  (c)  Zonal  velocity  u  and 
(d)  meridional  velocity  v  computed  on  NLOM  output  decimated  by  a  factor  of  8  in  both  directions  to  a  548 
by  1024  grid. 


grid  and  on  a  grid  that  is  decimated  by  a  factor  of  8  in 
both  directions,  to  a  548  by  1024  grid  that  is  comparable  in 
resolution  (1/4°  in  the  north-south  direction,  0.35°  in  the 
east-west  direction)  to  the  AVISO  grids.  In  the  case  of  the 
decimated  grid  only  the  subsampled  grid  points  are  used  in 
the  calculation,  so  that  the  grid  spacing  is  8  times  larger  than 
on  the  original  grid.  The  3-point  stencil  differences  are  sig¬ 
nificantly  different  from  7-point  and  9-point  differences,  and 
display  a  speed-dependent  bias,  on  the  high  resolution  grid  as 
well  as  on  the  decimated  grid.  At  some  grid  points  the  dif¬ 
ference  between  3 -point  and  9-point  estimates  are  as  large  as 
10%.  Results  from  the  stochastic  model  in  section  7  are 
consistent  with  this  suggestion  of  a  speed-dependent  bias  in 
3-point  stencil  results  even  in  the  limit  of  high  resolution.  The 
smaller  scatter  (i.e.  the  smaller  discrepancies  seen  between 
narrow  stencil  results  and  9-point  stencil  results)  seen  in 
Figure  4  relative  to  Figure  5  is  consistent  with  the  weak 
signal  in  the  AVISO  data  at  scales  smaller  than  the  2-3° 
wavelength  resolution  limitation  (Figure  1)  [Chelton  et  al.y 
2011]. 


[20]  In  Figure  6  we  display  the  zonally  averaged  dis¬ 
crepancies  between  the  squares  of  zonal  velocity  components 
computed  from  AVISO  using  stencils  of  differing  width: 


where  []  represents  a  zonal  average  operator.  We  also  display 


the  analogous  quantity  for  v  (meridional  velocity).  In  all 
cases  we  find  that  the  D  values  lie  much  closer  to  zero  for 
7-point  stencil  results  than  for  3-point  stencil  results,  with 
5-point  stencil  results  lying  in  between  these  extremes.  This 
demonstrates  a  “convergence”  in  the  calculation  of  velocities 
as  stencil  width  increases.  For  both  u  and  v  and  on  both  the 
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(a)  u,  1/3  degree  (b)  v,  1/3  degree 


3 

D 


0.1 

0 

-0.1 

-0.2 


-0.3 


-0.4 


-0.5 


(c)  u,  1/4  degree 


-50  0  50 

Latitude  (degrees) 


(d)  v,  1/4  degree 


Figure  6.  Values  of Du  and  D„  the  zonally  averaged  discre¬ 
pancies  between  the  squared  velocity  components  computed 
from  3-,  5-,  and  7-point  stencils  and  those  computed  from 
9-point  stencils,  (a)  Du  values  computed  on  the  1/3°  Mercator 
grid  reference  AVISO  product  for  November  5, 2008.  (b)  As 
in  Figure  6a  but  for  Dy.  (c)  As  in  Figure  6a  but  for  1/4°  latitude- 
longitude  grid  reference  AVISO  product  for  November  5, 
2008.  (d)  As  in  Figure  6c  but  for  D ^ 


1/3°  and  1/4°  grids,  the  3-  and  5-point  stencils  yield  velocity 
variances  that  are  too  weak  with  respect  to  the  9-point  sten¬ 
cils,  over  most  latitudes,  consistent  with  the  speed-dependent 
biases  seen  in  Figure  4.  On  the  isotropic  1/3°  Mercator 
grid,  the  values  of  Du  and  Dv  are  comparable,  with  typical 
values  of —0.1  for  the  3-point  stencil.  On  the  anisotropic  1/4° 
latitude-longitude  grid,  the  3-point  stencil  Du  values  are  as 
low  as  -0.2  to  -0.3  in  mid-  and  high-latitudes,  whereas  Dv 
values  seldom  drop  below  —0.1.  This  is  consistent  with  the 
fact  that  the  y-spacing  on  the  anisotropic  grid  remains  rela¬ 
tively  wide  at  high  latitudes,  unlike  the  x-spacing  at  high 
latitudes  on  the  anisotropic  grid,  or  the  spacing  in  either 
direction  at  high  latitudes  on  the  isotropic  grid  (recall  that  the 
error  estimates  increase  with  increasing  grid  spacing).  On  the 
anisotropic  1/4°  grid,  increasing  the  stencil  width  reduces  not 
only  the  D  values,  but  also  the  difference  between  values  of 
Du  versus  Dv. 

[21]  It  should  be  noted  that  some  of  the  differences 
between  results  on  the  1/4°  versus  1/3°  AVISO  grid  may  be 
due  to  the  bi-linear  interpolation  from  the  latter  to  the  former, 
as  well  as  to  the  differences  in  the  two  grid  spacings.  As 


noted  earlier,  the  bi-linear  interpolation  introduces  an  artifi¬ 
cial  sidelobe  into  the  zonal  wave  number  spectrum  of  kinetic 
energy  on  the  1/4°  grid.  To  investigate  the  effects  of  inter¬ 
polation  further,  we  have  interpolated  the  idealized  model 
snapshot  discussed  in  section  3  to  coarser  grids,  and  we  find 
that  the  difference  between  narrow  and  wide  stencil  estimates 
increases  over  that  found  on  the  original  higher  resolution 
grid.  Interpolation  of  the  idealized  model  output  to  a  coarser 
grid  which  anisotropic  (more  widely  spaced  in  one  direction 
than  the  other)  yields  a  greater  sensitivity  to  stencil  width  for 
derivatives  computed  in  the  widely  spaced  direction  than  for 
derivatives  computed  in  the  other  (better-resolved)  direction. 
All  of  this  is  consistent  with  the  results  in  Figure  6. 

[22]  Figures  7  and  8  display  Du  and  Dv  values  computed 
from  NLOM.  We  show  results  computed  on  the  original 
4384  by  8 1 92  grid  as  well  as  on  this  grid  decimated  by  factors 
of  2, 4,  and  8  in  both  horizontal  directions.  As  in  the  AVISO 
results,  the  difference  between  narrow  stencil  results  and 
9-point  stencil  results  decreases  rapidly  as  the  stencil  width 
increases,  for  both  velocity  components.  The  disparity 
between  3-point  stencil  results  and  wider  stencil  results 
increases  as  the  grid  spacing  increases.  The  results  for  both 


(a)  u,  native  NLOM  grid 


(b)  v,  native  NLOM  grid 


(c)  u,  NLOM  grid  decimated  by  2  (d)  v,  NLOM  grid  decimated  by  2 


Figure  7.  As  in  Figure  6,  but  using  February  15,  2002 
NLOM  output  rather  than  gridded  altimeter  data.  Values 
of  Du  and  D„  the  zonally  averaged  discrepancies  between 
the  squared  velocity  components  computed  from  3-,  5-,  and 
7-point  stencils  and  those  computed  from  9-point  stencils, 
(a)  Du  and  (b)  Dv  computed  from  NLOM  output  on  a  4384 
by  8192  grid,  (c)  Du  and  (d)  Dv  computed  on  NLOM  output 
decimated  by  a  factor  of  2  in  both  directions. 
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(a)  u,  NLOM  grid  decimated  by  4  (b)  v,  NLOM  grid  decimated  by  4 


(c)  u,  NLOM  grid  decimated  by  8  (d)  v,  NLOM  grid  decimated  by  8 


-50  0  50  -50  0  50 

Latitude  (degrees)  Latitude  (degrees) 


Figure  8.  As  in  Figure  7,  but  for  (a  and  b)  NLOM  grids 
decimated  by  a  factor  of  4  in  both  directions,  and  for  (c  and  d) 
decimation  by  a  factor  of  8  in  both  directions. 


Du  and  Dv  on  the  coarsest  grid  shown  (decimation  by  factor 
of  8;  Figures  8c  and  8d)  are  remarkably  similar  to  the  results 
on  the  AVISO  grids,  especially  the  AVISO  1/4°  grid.  In  all 
of  the  latter  cases,  the  3-point  stencil  Du  and  Dv  values  can 
reach  as  low  as  —0.1  to  —0.3.  Results  very  similar  to  those 
in  Figures  8c  and  8d  are  obtained  by  smoothing  the  high- 
resolution  NLOM  results  onto  a  grid  with  eight  times  lower 
resolution,  rather  than  by  subsampling  (decimation).  For  the 
sake  of  brevity  the  smoothed  results  are  not  shown  here. 
Having  used  the  high-resolution  model  to  demonstrate  that 
sensitivity  to  stencil  width  is  not  merely  an  artifact  of  the 
data  processing  inherent  in  the  creation  of  gridded  altimeter 
products,  we  focus  only  on  analysis  of  altimeter  data  in  the 
remainder  of  the  paper. 

5.  Effects  of  Non-uniform  Grid  Spacing 

[23]  The  formulae  in  section  2  have  assumed  that  the  grid 
spacing  as  measured  in  kilometers  is  uniform  along  a  par¬ 
ticular  direction  of  interest,  though  it  may  vary  between  the 
x-  and  y-directions.  However,  on  the  1/3°  AVISO  grid,  the 
grid  spacing  as  measured  in  kilometers  in  the  y-direction  is 
non-uniform.  Here  we  will  describe  a  method  which  accounts 
for  non-uniform  grid  spacing  such  as  this  in  computations 
using  classic  centered  differences. 

[24]  Suppose  we  are  interested  in  a  2m  +  1  point  stencil 
on  a  non-uniform  grid.  As  in  section  2,  let  the  central 
point,  where  the  derivative  is  to  be  estimated,  be  r.  Let  the 


coordinates  of  the  grid  points  in  the  stencil  be  r  +  hm, 
r  +  h„-  ,,  ...,r+hur,r-  h_{, r  -  A_(m  _  r  -  A_m. 
Let  us  expand  in  Taylor  series 


v(r  +  K)  =  r](r)  +  h.rf])(r)  + 


hW2)(r) 


2! 


+  ...+ 


K^\r) 


P- 


(16) 


and 


rfr  -  A_)  =  r Kr)  -  A_.f/'>(r)  +  hi^{r)  +  ... 

,  (~1  yV-rPir)  , 

P\ 

We  want  to  find  coefficients  cn  such  that 


(17) 


CmT){r+  hm)  +  cm-\Ti{r  +  hm-\)  +  ...  +  cxr](r  +  h\) 

+  c0r)(r)  +  c^rj{r  -  h-i)  +  ...  +  c_(m_!)r/(r  -  h_{m_ I}) 

+  c-mr](r-h-m)  ~r?(1)(r),  (18) 


to  as  high  an  order  as  the  stencil  width  allows.  For  example, 
for  a  3-point  stencil,  solving  the  matrix  equation 


yields  the  solution 

A-. 

C|  A,  (A, 

h\  —  A_| 

C0_  A.A-,  ’ 

_  -h\ 

"A_,(A,+A_,)’ 


(19) 


(20) 


(21) 


with  leading  order  error  term  These  formulae  col- 

lapse  to  the  uniform  grid  formulae  in  the  case  h\  =  A_|.  In  the 
y-direction  on  the  1/3°  AVISO  grid,  where  h{  and  /i_  j  are 
slightly  unequal,  c0  will  be  small  but  not  zero,  and  and  c_i 
will  be  nearly  equal  but  not  exactly  so.  For  a  5-point  stencil 
on  a  non-uniform  grid  the  coefficients  cn  are  given  by 


(  C2  \ 

1 

1 

1  1 

«  N 

(  o\ 

Cl 

hi 

A. 

0  —  A_i 

— A-2 

1 

Co 

C-l 

= 

A 

Wi 

0  Ai, 

0  -Al, 

hli 

-Aij 

0 

0 

V  C— 2  / 

A, 

0  Ai, 

A-2  ) 

w 

(22) 


The  7-  and  9-point  stencil  coefficients  on  a  non-uniform  grid 
are  solved  for  in  like  manner. 

[25]  We  have  computed  3-,  5-,  7-,  and  9-point  stencil 
derivatives  for  1/3°  AVISO  data  using  these  formulae  for  a 
non-uniform  grid.  The  results  are  not  shown  for  the  sake  of 
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(a) 

0.1 5 1 - 1 - - 


_  0.1 


-1  0  1 
£  from  9  point  stencil  (10*4  s'1) 


(b) 


Figure  9.  Effects  of  stencil  width  on  estimates  of  geostrophic  vorticity  £  computed  on  the  1/4°  latitude- 
longitude  grid  reference  AVISO  product  for  November  5, 2008.  (a)  Scatterplot  of  C  computed  from  3-point, 
5-point,  and  7-point  stencils  minus  £  computed  from  9-point  stencils  (y-axis),  plotted  versus  £  computed 
from  9-point  stencils  (x-axis).  Red,  green,  and  black  dots  denote  results  obtained  from  3-point,  5-point, 
and  7-point  stencils,  respectively,  (b)  Values  of  D ^  the  zonally  averaged  discrepancies  between  squared 
vorticities,  computed  from  3-,  5-,  and  7-point  stencils  and  those  computed  from  9-point  stencils. 


brevity.  The  differences  between  the  zonal  velocity  compo¬ 
nent  u  computed  using  the  non-uniform  grid  procedure  above 
versus  that  computed  from  uniform  grid  formulae  are  mea¬ 
surable  but  significantly  smaller  than  the  differences  arising 
from  different  stencil  widths.  The  equivalents  of  Du  and  Dv, 
for  example,  computed  from  differences  between  results 
using  the  non-uniform  grid  versus  uniform  grid  formulae,  are 
of  order  1%  for  the  AVISO  1/3°  Mercator  grid,  rather  than 
of  order  10-20%  as  seen  with  differing  stencil  widths.  This 
result  will  differ  for  more  rapidly  changing  grid  spacing.  For 
the  sake  of  simplicity  we  will  continue  to  use  the  uniform 
grid  spacing  formulae  in  the  remainder  of  the  paper. 

6.  Effect  of  Stencil  Width  on  Vorticity  Estimates 

[26]  In  this  section  we  briefly  discuss  the  effects  of  stencil 
width  on  estimates  of  vorticity  computed  using  classic 


centered  differences  from  the  1/4°  grid  AVISO  sea  surface 
height  reference  product  for  November  5,  2008.  In  Figure  9a 
we  show  a  scatterplot  of  C narrow  ~  C9 Pt  plotted  versus  (,9ph 
where  again  “narrow”  denotes  either  a  3-point,  5-point,  or  7- 
point  stencil  and  “9 pt”  denotes  a  9-point  stencil.  As  in  the 
scatterplots  of  velocity  (Figures  4  and  5),  the  7-point  stencil 
result  lies  closer  to  the  9-point  result  than  does  the  5-point 
stencil  result,  and  much  closer  than  does  the  3-point  stencil 
result.  In  Figure  9b  we  display 


where  Cnarrow  is  computed  on  the  narrower  3-,  5-,  and  7-point 
stencils,  and  £ 9pt  is  computed  on  the  9-point  stencil.  As  in 
the  plots  of  Du  and  Dv  (Figures  6-8),  we  see  a  decrease  to 
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near-zero  D  values  with  increasing  stencil  width.  Note  that 
omitting  the  (3  term  (the  last  term  on  the  right-hand  side  of 
(3))  in  the  estimate  of  C  yields  results  that  look  virtually 
unchanged  from  those  shown  in  Figure  9.  Thus,  as  expected, 
the  (3  term  does  not  play  a  first-order  role  in  the  computation 
of  vorticity. 

7.  A  Stochastic  Model  for  the  Speed-Dependent 
Bias  in  Derivative  Estimates 

[27]  In  Figures  4  and  5,  we  saw  that  the  difference  between 
“narrow”  stencil  and  9-point  stencil  estimates  of  geostrophic 
velocity  is  generally  positive  for  negative  values  of  u  or  v, 
and  negative  for  positive  values.  We  referred  to  this  in 
section  4  as  a  “speed-dependent  bias”.  Figures  6-8  demon¬ 
strate  that,  as  a  result,  values  of  Du  and  Dv  are  generally 
negative,  by  an  order  of  10%  for  a  stencil  width  of  three 
points.  Figure  9  demonstrates  that  the  same  principles  hold 
for  geostrophic  vorticity  estimates. 

[28]  Here  we  construct  a  stochastic  model  which  qualita¬ 
tively  captures  these  biases  in  the  narrow-stencil  estimates 
made  using  classic  centered  differences.  For  simplicity, 
we  discuss  here  only  the  first-derivative  (velocity)  results. 
We  also  consider  the  grid  spacing  h  to  be  constant  since  the 
results  of  section  5  show  that  the  errors  from  the  effects  of  the 
small  latitudinal  variation  of  grid  spacing  in  the  y-direction 
for  isotropic  grids  are  much  smaller  than  the  effects  arising 
from  the  coarseness  of  the  grid  spacing.  Our  stochastic  model 
is  a  simple  first  order  auto-regressive  model  (AR1).  We 
represent  77  as 


rj(r  4-  h)  =  <jyq{r)  -f  w(r  4*  h)}  (24) 

where  vv(r  -f  h)  is  a  value  taken  from  a  purely  random  process 
W,  also  called  a  “Gaussian  white-noise-in-space  stochastic 
process”.  That  is,  to  each  point  in  space  r  we  associate  a 
stochastic  process  W  giving  values  w(r)  that  are  completely 
uncorrelated  with  w(r  4  6)  at  a  neighboring  point  r  4  6  where 
6^0.  The  parameter  0  <  (j>  <  1  determines  the  spatial  auto¬ 
correlation  of  rj(r)  and  in  the  limit  <f>  — ►  0  we  see  that  77 
approaches  the  Gaussian  white-noise-in-space  stochastic  pro¬ 
cess.  A  similar  equation  to  (24)  applies  at  other  points.  In 
particular,  assuming  homogeneous  statistics,  we  can  write 


where 


1T=- +  --?*»  =^,(,1179^.  (29) 

\2h  12  h  h 

Here  we  have  used  the  fact  that  w(r  4  2 h)  and  vv(r  —  2 h)  are 
uncorrelated,  so  that  their  sum  is  another  random  function 
with  a  magnitude  obtained  from  quadrature  (i.e.,  taking  the 
square  root  of  the  sum  of  the  squares).  Inspection  of  this 
relationship  reveals  that  the  5-point  stencil  estimate  is  equal 
to  a  constant  ^2  times  the  3-point  stencil  estimate,  plus  a 
purely  random  process  W\  It  is  important  to  note  that  W'  is 
uncorrelated  with  the  3-point  stencil  estimate  since  it  is  the 
linear  combination  of  w  from  points  outside  the  3-point 
stencil.  Since  ^  exceeds  one,  the  N  =  3  derivative  displays 
a  speed-dependent  bias  relative  to  the  more  accurate  N  =  5 
stencil;  when  the  derivative  is  positive  the  discrepancy  is 
negative  and  vice  versa. 

[29]  Similar  manipulations  yield 


xfS  0)  9  r;(r  +  2h)  —  t?(r  —  2h) 

)  30%=3(  )  30  2h 

1  7?(r  +  3ft)  -  T?(r  -  3/i) 

30  2 h 


(30) 


and 


(„  672  (1)  1 68  jj(r  +  2h)  —  r)(r  —  2h) 

%=9 (n  ~  42q%=3W  -  420  Th 


4* 


32  r](r  4-  3 h)  —  rj(r—3h) 
420  Th 

3  T](r  4-  4 h)  —  77(7*  —  4 h) 
_  420  Th 

(672-  168<£4-3202-3<£3) 
420 


n£i  3(r)  +  w, 


(31) 


where 


rj(r  -  h)  =  077(r)  4-  w(r  -  h)  (25) 

r)(r  4-  2 h)  =  (jn](r  4-  h)  4*  w(r  4*  2h)  (26) 

r](r  —  2 h)  =  <fyq(r  —  h)  4-  w(r  —  2 h).  (27) 

Substituting  these  AR1  models  into  the  formulae  in  Table  1, 
we  find  the  following  relation  between  the  first  derivatives 
computed  with  a  3-point  stencil  (subscript  N  =  3  below)  and  a 
5-point  stencil  (subscript  N  =  5): 


W  —  — 
4* 


9  [w(r  +  2h)-w(r-2h)] 
60  h 

1  [w(r  4*  3 h)  —  w(r  —  3h)] 
60  h 


and 


8(1)  M  <t>  Hr  +  h)  -  v(r  -  h)] 
-  6%-3 W  -  6  2 h 

[w(r  4*  2h)  —  w(r  —  2h)\ 

12 h 


(28) 


Wm  =  - 

+ 


168  {w(r  +  2h)-w{r~2.h)\ 
840  h 

32  (w(r  +  3/i)-w(r-3A)] 
840  h 

3  [w(r  +  Ah)  —  w(r  —  Ah)) 
840  h 


7585)4  W 

-  —  —  =  0.2880—. 
840  h  h 


(32) 


(33) 
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Figure  10.  Normalized  difference  of  time-  and  zonally- 
averaged  zonal  and  meridional  velocity  variance,  [<  u2  -  v2  >]/ 
[<  u2  +  v2  >],  where  []  is  the  zonal  average  operator  and  <  > 
is  a  time-averaging  operator,  computed  from  16  years  of 
AVISO  products  (December  30,  1992  through  December 
31,  2008).  Both  the  isotropic  1/3°  Mercator  grid  reference 
AVISO  product  and  the  anisotropic  1/4°  latitude-longitude 
grid  reference  AVISO  product  are  used.  Values  of  time- 
averaged  u  and  v  are  removed  from  the  u  and  v  fields  before 
the  difference  is  computed.  Differences  are  computed  using 
(a)  3-point  stencils,  (b)  9-point  stencils,  and  (c)  a  blended 
product,  in  which  first  the  9-point  stencil  is  used,  then  a 
7-point  stencil  is  used  to  fill  in  missing  points  along  bound¬ 
aries,  followed  by  usage  of  a  5-point  stencil,  and  finally  a 
3-point  stencil. 


[30]  We  are  interested  in  the  speed-dependent  bias  in  the 
derivative  estimate  resulting  from  the  3-point  stencil  relative 
to  the  9-point  stencil,  so  we  write  (31)  in  the  form 

-  nJvisM  =  \  **"”>  (34) 

where  the  relative  bias  RB  =  *=£  and  c  is  the  coefficient  in 
front  of  7 fyL  3 (r)  in  (31): 

672  —  1 68</>  4-  32<f>2  —  3</>3 
C  = - 420 - '  (35) 


The  relative  bias  is  a  negative  and  increasing  function  on  the 
interval  0  ^  [0,  1).  The  bias  with  the  smallest  magnitude 
occurs  in  the  limit  of  very  strong  autocorrelation  0  — ►  1.  In 
this  limit,  we  have  (ignoring  the  random  parts) 

V/llji')  =  (36) 

(i)  /  x  533  (i) 

=  42o%i3M> 

which  yields,  for  the  differences  between  narrower  stencil 
derivatives  and  the  9-point  stencil  derivative  (again  ignoring 
the  random  parts), 

rJ/i)-s(r)  -  '&-<>('•)  =  (37) 

'JaJItW  -  V(N-9(r)  = 

[31]  It  is  important  to  note  that  the  random  parts  are  cor¬ 
related  with  the  discrepancy  between  the  3-point  and  9-point 
stencil  estimate  because  it  is  the  linear  combination  of  the 
white-noise  contribution  w  to  77  from  points  inside  the  9-point 
stencil.  Thus  we  do  not  expect  the  stochastic  model  to  give 
quantitative  predictions  of  the  speed-dependent  bias.  Indeed, 
the  stochastic  model  prediction  that  the  3-point  stencil  results 
will  be  biased  with  respect  to  the  9-point  stencils  by  about 
21%,  is  an  over-estimate.  Least-squares  fits  to  the  plots  in 
Figure  4  (based  on  altimetry  data)  yield  speed-dependent 
biases  in  the  range  4-7%,  and  least  squares  fits  to  the  plots  in 
Figures  5c  and  5d  (NLOM  decimated  by  a  factor  of  eight) 
yield  speed-dependent  biases  in  the  range  7-10%.  The  speed- 
dependent  bias  indicated  by  the  least  squares  fit  to  Figures  5a 
and  5b  (NLOM  on  its  native  high-resolution  grid)  is  much 
smaller  (1-2%).  Note  that  in  this  regard  the  appearance  of  the 
plots  in  Figure  5  is  deceptive.  There  are  many  more  points 
plotted  in  Figures  5a  and  5b  than  in  Figures  5c  and  5d.  Hence 
there  are  a  greater  number  of  points  with  large  misfits  in 
Figures  5a  and  5b;  however  the  percentage  of  points  having  a 
large  misfit  is  much  smaller.  The  smaller  bias  on  the  higher 
resolution  grid  is  consistent  with  the  discussion  of  the  von 
Neumann  analysis  in  section  3.  Clearly,  in  all  of  the  plots  in 
Figures  4  and  5,  discrepancies  between  narrow  stencil  esti¬ 
mates  and  9-point  stencil  estimates  at  individual  grid  points 
can  be  much  larger  than  the  least-squares  slopes  quoted 
above.  The  stochastic  model  appears  to  qualitatively  explain 
the  speed-dependent  bias  seen  earlier  in  the  paper,  though  the 
quantitative  estimate  of  the  bias  is  too  large.  The  stochastic 
model  indicates  that  this  speed-dependent  bias  will  be  seen 
even  if  derivatives  are  taken  at  very  small  grid  spacing.  This 
explains  why  a  speed-dependent  bias  was  seen  even  on  the 
native  high-resolution  NLOM  grid. 

8.  Effect  of  Stencil  Width  on  Estimates 
of  Anisotropy  in  Kinetic  Energy 

[32]  One  of  our  motivations  for  investigating  the  effects  of 
stencil  width  on  velocity  estimation  is  to  determine  whether 
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Table  2.  Formulae  for  Noise-Suppressing  Differentiators11 
A6  Formulae 


5 

7 

9 


7 

9 

5 

7 

9 


Low-Noise  Lanczos  Differentiator c 

2tKr+2* *)-hT?(r^/,)-^r-*)-2tKr-2A) 

USA 

lrf(r+'ih)+2Ti(r+2h)+i](r+h)-Ti(r-h)-2Ti(r-2h)-‘iT](r-‘ih) 

28  A 

4rKr-MA)-K3i7(r-t-3A)-t-2T?(r-t-2A)-|-fKr+*)-TKr-*)-2TKr-2A)-3rKr-3A)-4T7(r-4*) 


Super  Lanczos  Low-Noise  Differentiator c 

-22rKH-3A)+67Ty(r+2A)+58q(r+A)— 58q(r— A)— 67t^r— 2A)+22q(r~3A) 

-86fKr+4A)+142»7(r+3A)+l93»Kf+2A)+l26»7{r+A)- l26fKr-A)- 193»K''-2A)-I42r?{r -3A)+86t7(r-4A) 

u'8gA 


Smooth  Noise-Robust  (n  ~  2)  Differentiator c 

T(r+2A^27(^AC27(r-A)-7(r-2A) 

8A 


t?(r+3A)+47(r+2A)+57(r+A)-57(r. 


»Kr+4A)+6fKr+3A)+14fKr+2A)+14rKr+A)-14»Kr-A)-  14fKr-2A)-6rKr-3A)-rKr-4A) 

j'28A 


Smooth  Noise-Robust  (n  =  4)  Differentiator 

-5»Kr+3A)+12»j{'’+2A)+39T7(r+A)-39rKr~A)-l2rKr-2A)+5f7(r-3A) 

96A 

-2t;(r-»-4A)-q(r+3A)+16t7(r+2A)+27>7(r+A)-27r^;— A)-l6r^r-2A)+q{r-3A)+2f7(r-4A) 


"Source:  Holoborodko  (2011). 

*Yaluc  of  Nt  the  width  of  the  stencil. 
cFormulae  for  first  derivative. 


this  had  any  impact  on  the  estimation  of  the  anisotropy  of 
time-  and  zonally-averaged  surface  ocean  velocity  variance 
fields  by  Scott  et  al.  [2008].  In  this  section  especially  it  is  of 
interest  to  examine  results  computed  from  both  anisotropi- 
cally  and  isotropically  gridded  altimeter  products.  In  this  as 
in  earlier  sections  we  continue  to  focus  on  derivative  esti¬ 
mates  made  with  classic  centered  differences.  As  in  the  work 
of  Scott  et  al  [2008]  (see  their  Figure  9b),  we  compute  the 
quantity 


A/z2  = 


[<  u2  -  V2  >] 
[<  K2+v2>]’ 


(38) 


where  <  >  is  a  time-averaging  operator.  See  also  Ducet  et  al 
[2000]  and  Scharffenberg  and  Stammer  [2010],  who  per¬ 
formed  similar  computations  on  altimeter  data.  As  in  the 
work  of  Scott  et  al  [2008],  for  this  calculation  we  remove 
time  averages  from  u  and  v  before  computing  Mz 2.  Here  we 
compute  Mz2  over  16  years  of  satellite  altimetry  data,  from 
December  30,  1992  through  December  31,  2008.  Figure  10a 
reveals  that  3-point  stencil  computations  of  Mz2  yield  dif¬ 
ferent  results  in  mid-to-high  latitudes  on  the  1/4°  grid  than  on 
the  1/3°  grid,  despite  the  fact  that  the  1/4°  grid  AVISO 
product  is  derived  from  the  1/3°  product  by  simple  bi-linear 
interpolation.  The  discrepancy  between  results  on  the  two 
different  grids  is  expected  based  on  the  discussion  in  the 
previous  sections.  On  the  1/4°  grid,  the  bias  toward  smaller 
values  of  kinetic  energy  in  the  3-point  stencil  computation  is 
especially  pronounced  for  u  at  mid-to-high  latitudes  (com¬ 
pare  Figure  6c  to  Figures  6a,  6b,  and  6d).  At  these  latitudes, 
the  AVISO  regridding  from  the  original  Mercator  1/3°  iso¬ 
tropic  grid  to  the  1/4°  anisotropic  grid  apparently  degrades  u 
more  than  v.  Again,  this  is  because  u  is  based  on  derivatives 
in.y,  and  at  high  latitudes  the>>  spacing  is  coarser  on  the  1/4° 
anisotropic  grid  than  on  the  Mercator  1/3°  isotropic  grid. 


Since  the  Jt-spacing  is  always  less  on  the  1/4°  anisotropic  grid 
than  on  the  Mercator  1/3°  isotropic  grid,  v  does  not  suffer 
from  this  degradation. 

[33]  Consistent  with  earlier  discussions,  Figure  10b 
demonstrates  that  Mz2  computed  with  9-point  stencils  yields 
very  similar  results  when  computed  on  the  1/4°  grid  as  on  the 
1/3°  grid.  However,  a  disadvantage  of  utilizing  wider  stencils 
is  that  more  information  is  lost  along  the  boundaries  where 
a  full  span  of  the  stencil  does  not  exist  because  of  land.  To 
alleviate  this  problem,  we  have  also  created  a  blended 
velocity  estimate,  which  begins  with  velocities  computed 
using  9-point  stencils,  then  where  possible  fills  in  missing 
velocities  computed  using  7-point  stencils,  followed  by  those 
computed  from  5-point  stencils,  and  finally  by  those  com¬ 
puted  from  3-point  stencils.  This  blended  estimate  has  the 
disadvantage  that  the  “quality”  of  derivatives  is  not  uniform 
across  all  grid  points,  and  the  advantage  that  the  number  of 
grid  points  for  which  a  velocity  estimate  is  available  is  as 
large  as  in  the  3-point  estimates.  In  Figure  10c  we  show  Mz2 
computed  from  the  blended  estimate.  As  in  Figure  10b,  the 
differences  between  Mz2  values  computed  on  the  1/3°  versus 
1/4°  AVISO  grids  are  much  smaller  than  seen  in  Figure  10a 
(3-point  stencil  computation).  Discrepancies  between  1/3° 
and  1/4°  results,  and  between  pure  9-point  stencil  results  and 
results  from  the  blended  product,  are  largest  at  high  southern 
and  northern  latitudes. 

[34]  Scott  et  al  [2008]  and  Scharffenberg  and  Stammer 
[2010]  both  found  near  isotropy  in  the  mid-  and  high- 
latitude  zonally  averaged  values  of  the  difference  between 
zonal  and  meridional  velocity  variances.  Scott  et  al  [2008] 
utilized  the  Mercator  1/3°  isotropic  grid  AVISO  product 
while  Scharffenberg  and  Stammer  [20 1 0]  utilized  data  from 
the  tandem  phase  of  the  JASON/TOPEX  missions.  Figure  10 
shows  that  near-isotropy  (small  values  of  Mz2)  is  also  seen 
when  wider  stencils  are  used.  We  thus  conclude  that  the 
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Figure  1 1 .  Response  of  an  ideal  differentiator  and  of  5-,  7-, 
and  9-point  stencils,  as  a  function  of  wave  number,  on  1/3° 
AVISO  grid,  for  (a)  Low-noise  Lanczos  (LNL)  differentiator 
and  (b)  smooth  noise-robust  (SNR;  n  =  2)  differentiator.  The 
grid  spacings  used  are  37  km,  and  the  response  is  the  absolute 
value  \A\  of  the  eigenvalue  of  the  von  Neumann  analysis. 
See  text  for  more  details.  Vertical  lines  are  drawn  at  the  wave 
number  (40  km)" 1 ,  corresponding  approximately  to  the  2-3° 
wavelength  resolution  limitation  of  the  AVISO  sea  surface 
height  fields. 


inference  of  Scott  et  al.  [2008]  and  Scharffenberg  and 
Stammer  [2010],  that  the  zonally  averaged  oceanic  kinetic 
energy  is  nearly  isotropic  at  midlatitudes,  is  not  an  artifact  of 
the  three-point  stencil  used  to  compute  derivatives  in  these 
studies. 

9.  Effects  of  Noise 

[35]  As  stated  by  Holoborodko  (2011),  in  many  applica¬ 
tions,  ideal  differentiators,  and  classic  centered  differences 
which  closely  approximate  ideal  differentiators  over  a  wide 
range  of  scales,  are  problematic.  This  is  because  in  many 
applications  there  is  significant  noise  at  high  frequencies 
(or  wave  numbers),  and  the  high  frequency  noise  should  be 
suppressed  if  one  wants  reasonable  estimates  of  the  deriva¬ 
tive  at  low  frequencies.  A  common  procedure  for  doing  this 
is  to  first  smooth  the  data  with  a  least-squares  polynomial  fit 
and  to  then  estimate  the  derivative  using  this  polynomial.  An 
extensive  literature  on  smoothing  and  differentiation  through 
least-squares  fitting  exists  in  the  chemical  and  signal  pro¬ 
cessing  community  [Savitzky  and  Golay ,  1964;  Steiner  et  al., 
1972;  Gorry ,  1990;  Luo  et  al.y  2005].  We  found  the  discus¬ 
sion  of  Holoborodko  (201 1)  especially  useful. 

[36]  Following  the  discussion  of  Holoborodko  (2011), 
we  considered  a  particular  class  of  Savitzky-Golay  filters, 
known  as  “low-noise  Lanczos”  and  “super  Lanczos  low- 
noise”  differentiators,  as  well  as  a  class  of  filter  designed  by 
Holoborodko  (201 1)  and  known  as  a  “smooth  noise-robust” 
differentiator.  Table  2  lists  the  first-derivative  formulae  for 
some  of  these  filters.  As  a  check,  we  derived  these  formulae 
for  ourselves,  and  obtained  the  same  answers  as  Holoborodko 
(2011). 

[37]  In  Figure  1 1,  we  display  results  of  the  von  Neumann 
analysis  for  the  low-noise  Lanczos  and  smooth  noise-robust 
(n  =  2)  differentiators.  See  Holoborodko  (2011)  for  a  full 
exposition  and  derivation,  which  includes  an  explanation  of 
the  meaning  of  n  in  the  smooth  noise-robust  differentiators. 
As  shown  in  the  figure,  and  as  described  by  Holoborodko 
(2011),  these  differentiators  are  designed  to  suppress  noise 
at  high  wave  numbers,  while  remaining  close  to  the  ideal- 
differentiator  behavior  at  low  wave  numbers.  However,  a 
problem  for  these  low-noise  differentiators  is  that  they  devi¬ 
ate  significantly  from  the  ideal  differentiator  response  over 
about  half  of  the  range  of  wave  numbers  lower  than  the  2-3° 
wavelength  resolution  limitation  of  the  AVISO  sea  surface 
height  fields.  Thus,  these  low-noise  differentiators  suppress 
much  of  the  signal  as  well  as  the  noise.  The  spectra  shown 
in  Figure  1  demonstrate  that  there  is  little  variance  in  the 
AVISO  data  in  wave  numbers  higher  than  those  representing 
the  resolution  limitation.  Thus,  for  the  applications  relevant 
to  this  paper,  noise  suppression  at  high  wave  numbers  may 
be  less  important  than  retaining  near-ideal  differentiation  at 
low  wave  numbers. 

[38]  To  test  this  further,  we  looked  at  differences  between 
5-point,  7-point,  and  9-point  stencil  results  with  these  low- 
noise  differentiators  applied  to  a  snapshot  of  the  1/3°  AVISO 
data  set  (Figure  12).  The  discrepancies  between  derivatives 
computed  with  narrow  and  wide  stencils  is  much  greater 
using  these  low-noise  differentiators  than  when  using  classic 
centered  differences.  We  obtained  similar,  though  some¬ 
what  less  extreme,  results  using  the  “super  Lanczos  low- 
noise”  differentiators  and  the  “smooth-noise  robust  (n  =  4)” 
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(a)  Centered  differences 
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Figure  12.  Difference  between  zonal  velocity  u  computed  with  5-point  and  9-point  stencils  (green  dots), 
and  with  7-point  and  9-point  stencil  results  (black  dots),  for  (a)  classic  centered  difference  formulae  used 
elsewhere  in  the  text,  (b)  low-noise  Lanczos  differentiators,  and  (c)  smooth  noise-robust  (n  =  2)  differentia¬ 
tors.  Differences  are  computed  for  the  November  5,  2008  snapshot  of  the  1/3°  AVISO  data  set. 


(c)  Smooth  noise-robust 
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differentiators  (results  not  shown).  If  we  accept  “conver¬ 
gence”  with  increasing  stencil  width  as  an  important  crite¬ 
rion,  we  conclude  that  the  low-noise  differentiators  are  not  as 
well-suited  as  classic  centered  differences  for  computing 
derivatives  on  current-generation  altimeter  data.  For  next 
generation  wide-swath  satellite  altimeter  data  [Fu  and 
Ferrari ,  2008],  the  noise  characteristics  of  the  signals  are  as 
yet  unknown.  Although  the  instrument  is  expected  to  per¬ 
form  some  amount  of  averaging  over  the  raw  measurements, 
the  amount  of  smoothing  will  be  much  less  than  that  used  in 
the  construction  of  current-generation  gridded  satellite 
altimeter  data.  As  a  result,  noise-suppressing  differentiators 
are  likely  to  become  more  necessary  and  important  with  the 
coming  of  high  resolution  wide-swath  data. 

[39]  We  also  investigated  the  extent  to  which  derivative 
estimates  made  with  classic  centered  differences  are  affected 
by  noise  in  the  altimeter  data.  For  this  purpose  we  have 
constructed  three  types  of  random  noise  fields  on  the  1/3° 
AVISO  grid.  The  first  is  simply  a  noise  field  which  is  random 
from  one  point  to  another.  The  second  and  third  versions  are 


produced  from  Blackman  filtering  applied  to  the  random 
noise.  The  second  version  utilizes  a  Blackman  filter  which 
goes  to  zero  at  a  distance  of  2°,  while  the  third  goes  to  zero  at 
3°.  Each  of  the  three  fields  is  multiplied  by  constants,  such 
that  we  end  up  with  fields  having  RMS  values  of  0. 1 , 0.5,  and 
1  cm.  This  yields  9  random  noise  fields  in  all.  In  Figure  13  we 
display  differences  in  zonal  velocity  estimated  with  9  point 
stencils,  with  and  without  the  RMS  1  cm  amplitude,  3°  fil¬ 
tered  noise.  We  compare  these  differences  (black  dots)  to  the 
differences  between  zonal  velocity  estimated  using  3  point 
versus  9  point  stencils  (red  dots,  no  extra  noise  included). 
The  noise  adds  an  uncertainty  in  the  velocity  estimates  which 
is  comparable  in  magnitude  to  the  differences  one  sees  in 
estimates  made  with  narrower  stencils.  Obviously,  the  size  of 
the  noise-induced  error  decreases  with  decreasing  RMS 
amplitude.  The  error  increases  with  decreasing  horizontal 
scale;  the  error  with  2°  noise  is  larger  than  with  3°  noise,  and 
the  error  with  random  spatial  noise  is  larger  still  (results  not 
shown).  However,  in  contrast  to  the  speed-dependent  bias 
one  sees  with  narrower  stencils,  the  uncertainty  added  by 
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Figure  13.  Difference  between  zonal  velocity  u  computed  from  the  November  5,  2008  snapshot  of  the 
1/3°  AVISO  data  set  with  3-point  minus  9-point  stencils  (red  dots),  and  between  u  computed  with  9-point 
stencils  with  1  cm  RMS  random  noise  filtered  to  3°  added  to  AVISO,  minus  u  computed  without  the 
added  noise  (black  dots). 


these  noise  fields  is  not  biased  in  an  obvious  way.  Instead  it 
appears  to  be  distributed  evenly  around  zero,  with  larger 
excursions  occurring  at  smaller  velocities.  In  addition,  the 
effects  of  the  noise  are  just  as  large  with  3-point  stencils  as 
with  9-point  stencils  (not  shown).  Therefore  the  wider  sten¬ 
cils  are  not  a  disadvantage  in  this  exercise.  The  differences  in 
estimates  of  the  derivatives  of  the  noise  fields  decrease  with 
increasing  stencil  width,  just  as  differences  in  the  derivative 
estimates  of  the  AVISO  signals  do  (not  shown). 

10.  Summary  and  Discussion 

[40]  Drawing  on  the  well-developed  numerical  methods 
literature,  we  have  shown  that  classic  centered-difference 
estimates  of  geostrophic  velocities  and  vorticities  from  grid- 
ded  products  of  satellite  altimeter-derived  sea  surface  height 
lie  closer  to  results  obtained  with  an  “ideal  differentiator”  as 
stencil  width  increases.  This  is  especially  true  for  mid-to-high 
latitude  velocities  computed  on  the  1/4°  AVISO  grid,  which 
is  anisotropic  (has  different  grid  spacings  in  the  east-west 
versus  north-south  directions).  Similar  computations  using 
NLOM,  a  high-resolution  numerical  model,  demonstrate  that 
the  impact  of  stencil  width  on  geostrophic  velocity  estimation 
is  not  an  artifact  of  the  particular  procedures  and  errors 
involved  in  creating  gridded  satellite  altimeter  products.  A 
stochastic  model  developed  here  qualitatively  explains  the 
speed-dependent  biases  seen  in  velocity  estimates  made  with 
narrow  (3-point)  stencils. 

[41  ]  This  study  was  inspired  in  part  by  our  earlier  study  of 
the  anisotropy  of  ocean  surface  velocity  variances  [Scott 
et  al. ,  2008],  which  utilized  the  traditional  thin  (3-point) 
stencils  on  the  1/3°  Mercator  AVISO  grid.  We  have  shown 
here  that  usage  of  the  3 -point  stencils  which  are  widely  used 
in  the  oceanographic  literature  gives  incorrect  estimates  of 


the  anisotropy  of  velocity  variance  on  the  1/4°  anisotropic 
AVISO  grid;  we  recommend  the  use  of  7-point  stencils 
instead.  Although  the  analyses  in  sections  2-7  show  clearly 
that  applications  that  require  precise  geostrophic  velocity 
estimates  should  utilize  wider  stencils,  the  results  of  section 
8  show  that  the  general  conclusions  of  Scott  et  al.  [2008] 
and  Scharffenberg  and  Stammer  [2010]  that  zonally  aver¬ 
aged  velocity  variance  in  midlatitudes  is  nearly  isotropic  still 
pertain  when  the  derivatives  arc  computed  more  accurately 
with  wider  stencils.  Our  study  also  has  relevance  to  the  pro¬ 
posed  future  wide-swath  satellite  altimeter  mission  [Fu  and 
Ferrari ,  2008],  which  is  expected  to  map  sea  surface 
heights  at  resolutions  about  20-30  times  higher  than  the 
AVISO  sea  surface  height  fields.  Consistent  with  predictions 
from  our  stochastic  model,  our  analysis  of  high  resolution 
NLOM  output  indicates  that  even  at  high  horizontal  resolu¬ 
tion  (small  grid  spacing)  there  are  significant  differences 
at  individual  grid  points  between  velocities  computed  via 
3-point  stencils  and  wider  stencils.  Thus  wider  stencils 
should  be  of  interest  for  future  satellite  altimeter  missions  as 
well  as  for  present  ones.  In  future  altimeter  missions,  the  data 
may  begin  to  approach  what  is  the  norm  in  many  other  fields, 
that  is,  a  well  resolved  low-wavelength  signal  with  sig¬ 
nificant  noise  at  higher  wavelengths.  In  this  case,  noise¬ 
suppressing  differentiators  such  as  those  used  extensively  in 
the  chemical  and  signal  processing  communities  will  become 
of  greater  interest  to  oceanographers. 
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